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Abstract 

A class of multivariate mixed survival models for continuous and 
discrete time with a complex covariance structure is introduced in a 
context of quantitative genetic applications. The methods introduced 
can be used in many applications in quantitative genetics although the 
discussion presented concentrates on longevity studies. The frame- 
work presented allows to combine models based on continuous time 
with models based on discrete time in a joint analysis. The continuous 
time models are approximations of the frailty model in which the haz- 
ard function will be assumed to be piece-wise constant. The discrete 
time models used are multivariate variants of the discrete relative risk 
models. These models allow for regular parametric likelihood-based 
inference by exploring a coincidence of their likelihood functions and 
the likelihood functions of suitably defined multivariate generalized 
linear mixed models. The models include a dispersion parameter, 
which is essential for obtaining a decomposition of the variance of the 
trait of interest as a sum of parcels representing the additive genetic 
effects, environmental effects and unspecified sources of variability; as 
required in quantitative genetic applications. The methods presented 
are implemented in such a way that large and complex quantitative 
genetic data can be analyzed. 
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1 Introduction 



Longevity is an important trait often considered in animal breeding programs 
[31, 30, 33, 37, 27, 35, 17, 9, 6, 32, 1]. Even small changes in the longevity 
of a population under production might have remarkable economic, welfare 
and ethics consequences [30, 12]. Since the study of longevity involves several 
types of incomplete observation (e.g. censoring, truncation, late entry and 
competing risks), survival and event-history- analysis techniques are typically 
used [8, 19, 3]. However, the use of those techniques in the context of quan- 
titative genetics of longevity involves several non-trivial challenges. We will 
present a model framework here that allows to overcome these challenges. 
These models will have a structure of means and covariances similar to the 
gaussian linear mixed models classically used in quantitative genetics. They 
will allow for a proper representation of quantitative genetic phenomena and 
for efficient implementations required in practice. We will show that these 
models extend the class of models currently in use for studying the genetic 
of longevity. 

The first challenge in the use of survival models for characterizing the 
genetic aspects of longevity is the high complexity of the models. Typically 
it is necessary to adjust simultaneously for the effects of several explana- 
tory variables, some continuous with linear effects (e.g. breed composition 
and heterosis [21]), some factors with many levels (e.g. herd) or even time- 
dependent effects (e.g. year and season). Further complexity is added by the 
necessity of representing complex genetic effects (several generations deep 
pedigrees) under different genetic models (e.g. sire model, sire-dam- model 
etc). A typical scenario of applications in quantitative genetics involves a 
very large number of observations (individuals); indeed often the analyses 
involve several hundred thousand of observations (we will present two rela- 
tive small illustrative examples involving 142,133 and 200,084 individuals). 
This leads to inference problems of the complexity equivalent to solve sys- 
tems of linear equations with several hundred of thousand, or even several 
million, simultaneous linear equations (in our examples between 8,248 and 
118,432 simultaneous linear equations) when using classic gaussian linear 
mixed models. Now, the complexity of survival and event-history models is 
even larger since these models require to keep track of the individuals' his- 
tory. Moreover, as it will be apparent from the description of the methods 
available and from our discussion, some of those models might present very 
flat likelihood functions making the inference problem hard and numerically 
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unstable. This rules out the naive use of the standard methods of survival 
analysis and makes the use of specially constructed models and approxima- 
tions mandatory. 

A second challenge is the proper representation of the genetic phenomena 
in play. The interpretation of the gaussian linear mixed models in genetic 
terms requires a linear decomposition of the so called total variance (i.e. the 
phenotypic variance after having corrected for the effects of a range of non- 
genetic related aspects with known effects on the trait of interest) [18, 34]. 
This decomposition of the total variance occurs naturally in the gaussian 
linear mixed models but not in all the adaptations of survival models spe- 
cially constructed for genetic evaluations of longevity. We will show that a 
way to circumvent this problem is to consider models containing a dispersion 
parameter which plays the role of the residual variance (i.e. the part of the 
total variance that is neither attributed to genetic components nor to iden- 
tifiable sources of variation) and base the inference on the quasi-likelihood 
theory [36, 4]. The use of generalised linear mixed models with a dispersion 
parameter occurs only occasionally in the literature of quantitative genetics 
[14, 5], and the genetics and the statistics consequences of using a free disper- 
sion parameter has not being systematically explored yet. Moreover, survival 
models with dispersion parameter were, to the best of our knowledge, never 
considered in the literature of quantitative genetics. 

A third challenge is the presence of incomplete observations following 
complex patterns. Variants of the so called Cox proportional model [8] were 
used for univariate traits describing longevity [9, 13, 10]. These techniques 
succeeded to implement models of approximations that were operational for 
some of the purposes of animal breeding. However, all those models were 
of univariate nature (i.e. they consider one trait at time) [30, 11], which 
makes it difficult, if not impossible, to address some key questions related 
to the issues of competing risks, informative censoring and inhomogeneity 
of the failure mechanism. A clear example is the study of time to death 
were the individuals on the study may die from one among several different 
causes of death. Moreover, the use of multivariate models is necessary to 
study the time death in a competing risk scenario, by evaluating the risk of 
death for each specific cause simultaneously. This allows to form the basis to 
study further aspects as the presence of common genetic determining factors 
and/or independent or specific genetic factors for each of the death causes. 
In the framework presented here, it will be possible to treat these issues by 
using suitable multivariate models. 
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Two types of models will be considered: models based on continuous time 
and models based on discrete time. The continuous time models (CTM) will 
be suitable approximations of the semi-parametric Cox proportional model in 
which the hazard function will be assumed to be piece-wise constant. These 
models allow for regular parametric likelihood-based inference by exploring 
a coincidence of their likelihood functions and the likelihood function of a 
Poisson model applied to a specially constructed data [2]. When including 
gaussian random components (log gaussian frailties) the inference can be 
based on a Laplace approximation to a high dimensional integral [4, 24, 28] . 
Here, we will include also a dispersion parameter that will play the role of the 
residual variance which is an essential element for the genetic interpretation 
of the models. The CTM will depend on an arbitrary choice of time cut 
points used to define the intervals at which the hazard function is assumed 
to be constant. 

The discrete time models (DTM) used are multivariate variants of the 
classic discrete relative risk models [19] . These models present advantages in 
terms of computational and statistical complexity as compared to the CTM: 
the algorithms run faster relatively to the analogous algorithm for CTMs, 
are numerically more stable, and the statistical inference is more efficient. 
The loss of genetic information occurred when switching from the CTM to 
the DTM is minimal in typical applications in quantitative genetics and this 
loss is fairly compensated by the avoidance of capturing a large amount of 
noise. Here, we will also introduce the use of a dispersion parameter, which 
will be essential for well representing genetic scenarios. 

The aim of this article is to introduce, characterize and discuss a class of 
multivariate mixed survival models for continuous and discrete time with a 
complex covariance structure in a context of quantitative genetics applica- 
tions. The methods introduced here can be used in many other applications 
in quantitative genetics although the discussion presented concentrates on 
longevity. 

The paper is organized as follows. Section 2 describes the basic set-up and 
the genetic scenario discussed. There, a suitable multivariate version of the 
proportional hazard model is introduced in general terms. Those models will 
encompass models for competing risks possibly defined with different types 
of time scale (continuous and discrete time). The techniques for statistical 
inference under those models are presented in section 3, including a general 
discussion on the calculations involved in the likelihood based inference and 
some connections with multivariate generalized linear mixed models (section 
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3.1). A Poisson approximation useful for efficiently implementing likelihood 
based inference for continuous time models is presented in section 3.2. Sec- 
tion 3.3 presents an extension of the models to a situation where there is a 
stratification variable. A discussion comparing the amount of statistical in- 
formation between models constructed with discrete time and models based 
on continuous time is given in section 3.4. The quantitative genetic theory 
behind the models considered here requires a special decomposition of the 
phenotypic variance in terms of the variance of the random components and 
a scale parameter present in the models. This will be necessary for the cal- 
culation of the so called heritability, which is crucial in quantitative genetic 
applications. Section 4 will discuss methods for that and some technical de- 
tails involving counting processes and martingale theory are presented in the 
appendix. Two illustrative examples involving longevity of sows and dairy 
cattle will be presented in section 5.1 and 5.2 respectively. Some discussion 
will be given in section 6. 

2 The basic set-up and genetic scenario 

The class of models we will discuss are thought to be applied in the following 
general scenario. The life spans of a population of n individuals are observed. 
Genetic information (typically in the form of a reasonably deep pedigree or 
comprehensive genomic data) and the values of k explanatory variables are 
available for each individual. The interest lies in modeling the longevity, 
i.e. the length of the life span or the length of the productive life of the 
individuals, with particular interest on some forms of genetic determination 
typically used in quantitative genetics. 

The longevity can be operationally measured using a continuous time, 
by determining the time elapsed between two life events in the life of the 
individuals {e.g. the time elapsed from the first parity the to the culling of 
cows under production) or a discrete time, by counting a number of certain 
events during the life of the individual {e.g. by counting the number of 
survived parities of cows). We will present two illustrative examples: one 
involving the longevity of sows and the other studying the longevity of dairy 
cattle. Many similar examples can be found in the literature with other 
husbandry animals as ewes [25], salmon [20] among others. 

An additional issue that might occur in longevity studies is the presence 
of competing risks, which arises when the individuals could be culled for 
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one of two or more types of culling. For instance, in the example of dairy 
cattle longevity, the cows might die or be slaughtered, and the interest is in 
studying the culling rates for both causes. We assume in this scenario that 
there are two causes of death or culling, which will be labeled by the index j. 
It is straightforward to extend the setup described here to the general case 
with more than two culling causes or to the case with only one culling type. 

The longevity will be characterized by the time development of the rate 
at which the individuals die or are culled. This is measured differently when 
using continuous or discrete time as we define below. In order to describe the 
models we have in mind, consider two random variables: T representing the 
time (continuous or discrete) of culling and J indicating the cause (or type) 
of culling. In the case of the continuous time, we define the cause- specific 
hazard function [19] for the j th culling cause (for j = 1, 2) by 

A,| W= li m P ^ r<t+ A A ' J = ^ r>t > ,V«>0. 

In contrast, when the time is discrete, we define the cause-specific hazard 
probability function [19] for the j th culling cause (for j = 1, 2) by 

\ m (t) = P(T = t,J = j\T >t) fort = 0,1, 2,... . 

Although these two characterisations of the time development of the culling 
rate are of different nature, they will essentially play the same rule in the 
models discussed and will be generically referred simply as the cause-specific 
hazard function. Note that in the case where only one cause of culling is stud- 
ied the cause-specific hazard function and the cause-specific hazard probabil- 
ity function defined above coincide with the hazard function and the hazard 
probability function used in the literature of survival analysis (see [19]). 

The models presented below will assume a scenario where there are two 
causes of culling (indexed by j). Here, the cause-specific hazard functions (or 
the hazard probability functions) are specified in terms of k explanatory vari- 
ables (the fixed effects) and a set of gaussian random components. Without 
loss of generality, we consider only two random components for each culling 
cause: Ui and U2 that will represent additive genetic effects (usually deter- 
mined using information on the pedigree of the animals in play) for cause 1 
and 2 respectively and Vi and V 2 , representing an environment effect {e.g. 
the effect of herd where an animal is kept) for cause 1 and 2, respectively. 
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The random components related to genetics (i.e. Ui and U 2 ) are indepen- 
dent of the components related to environment (Vi and V2) by construction. 
This constraint can be relaxed for describing possible interactions between 
environmental and genetic effects, but we will not pursue this project here. 
The extension of the present definition to a situation with a different number 
of random components is straightforward. According to the general model 

t h 

we describe, the vector of cause specific hazard functions of the i individual 
(for % — 1, . . . ,n), conditionally on Ui = Ui, U2 = u 2 , Vi = Vi and V 2 = v 2 
(for the realizations u x , u 2 , Vi and v 2 of Ui , U 2 , V 2 and V 2 , respectively), 
is given by 



Here the times t\ and t 2 might be continuous or discrete; moreover, t\ and 
t 2 are not necessarily of the same type (as in the example of sows longevity 
considered in section 5.1). Equation (1) holds for tj > when tj is continuous 
and for tj = 0,1,2,... if tj is discrete (j = 1,2). Furthermore, Ai(-) and 
A 2 (-) are baseline hazard functions describing a common time development 
of the cause specific hazard functions (typically considered as a nuisance 
parameter), (3 1 ,(3 2 £ K fc are (finite dimensional) parameters representing 
the fixed effects (also viewed as a nuisance parameter), Xj, Zji and Wjj are 
incidence matrices for the fixed, genetic and environment effects for the 
cause and the i th individual. When the time is continuous \{j],i(-) in (1) is 
a hazard function and the marginal model described will correspond to a 
variant of the Cox proportional hazard model with frailties [19, 3]. If the 
time is discrete, A[i] ;i in (1) is a hazard probability function and the marginal 
model described will be a variant of the discrete relative risk model with 
frailties [19]. 

In the case of the models based on continuous time, we assume that the 
baseline function, Xj(-) is piece-wise constant. That is, we assume that there 
is a set of disjoint intervals, I±, . . . , Ik covering the positive real line, in which 
Xj(-) is constant. This constraint in the class of possible baseline functions 
will allow us to connect the models described here with some generalized lin- 
ear mixed models: Moreover, this will allow us to perform efficient inference 
with complex data using existent software. 

The structure of covariance of the random components is a crucial part 



A[i ])i (ti|Ui,Vi) 
A [2] ,^2|U 2 ,V 2 ) 
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of the construction of the models for longevity described since this is the 
part of the model where the additive genetics and the environment con- 
tributions to the cause specific hazard functions are considered. We as- 
sume the random components to be multivariate normally distributed, i.e. 
(Ui, U 2 , Vi, V 2 )' ~A/"(0,£) where 



£ = 



Si <g> A 
£ 2 (8)1 



(2) 



Here, I is an identity matrix and A is the n x n additive relationship matrix 
constructed as follows. Each entry of A represents a pair of individuals in 
the population. An — 1 for i — 1, . . . , n, and Aij = if the pair of individuals 
(i, j) are not related. In the case were the pair of individuals (i, j) are related, 
we consider the degree of relationship r given by 1 for relatives of first-degree, 
2 for relatives of second- degree, 3 for relatives of third-degree and so on. In 
that case = (l/2) r . The additive relationship between two individuals 
measures the proportion of identical by descent genes expected to be shared 
by pairs of individuals in the population. The methods introduced here 
could also be applied to genomic analyses where the relationship matrices 
are inferred from genomic data. 

The covariance structure between the genetic random components is given 

by 



Tg.12 



(Jg.11 



In particular the genetic correlation between the two culling causes is given 



by Pg = 



random components is given by 



Finally, the covariance structure between the environment 



u e.l 

c e .i2 cr 



Oe.12 
2 
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3 Inference 



3.1 Representation of the counting process and con- 
nections with generalized linear mixed models 

The statistical inference for the longevity models can be made by using the 
coincidence of the conditional likelihood of those models (conditioned on the 
frailties) with the conditional likelihood of a generalized linear model applied 
to a specially designed pseudo data representing an underlying counting pro- 
cess as briefly described below. In the case of the continuous time, the pseudo 
data is constructed by generating one pseudo observation for each time inter- 
val (in which the baseline function is assumed to be constant) observed for 
each individual. In the case of the discrete time, the pseudo data contains one 
observation per each unit of time observed for each individual. To facilitate 
the description we use the term "observed times" to refer to the observed 
discrete times or the time intervals related to the observed continuous time. 
A pseudo response variable was created for each of the culling causes. These 
variables are such that they take the value when an individual is at risk of 
being culled for the current reason but was not culled or take the value 1 if 
the individual was culled for the current reason in that period. This pseudo 
data describes the groups of individuals at risk at each observed time. 

The statistical inference is performed by using a multivariate generalized 
linear mixed model (MGLMM) applied to the pseudo data. One dimensional 
versions of the calculations we describe below were given previously in [2] for 
models without random components and in [15] for mixed models for con- 
tinuous time models (piece-wise constant baselines). For one-dimensional 
discrete times models without random components see [19]. Each of the 
marginal model of this MGLMM corresponds to a culling reason and is de- 
fined with a logarithmic link function. The distribution of each marginal 
model is Bernoulli (or binomial) when the corresponding time is discrete or 
Poisson when the time is continuous. Furthermore, the marginal models con- 
structed with continuous time include an offset variable with the logarithm 
of the length of the corresponding time intervals. Additionally, a discrete ex- 
planatory variable counting the order of the time period for each individual 
should be included in the model in order to represent the baseline function. 
The model should also specify the covariance structure of the random com- 
ponents given in equation (2). Furthermore, each marginal model includes a 
dispersion parameter <pj (known in the literature of generalized linear mod- 
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els as an over-dispersion parameter for binomial and Poisson models), which 
allowed us to better characterizing the genetic scenario. 

The MGLMM described above can be adjusted with any software im- 
plementing multivariate generalized linear mixed models, which allows the 
inclusion of over-dispersion parameters in binomial and Poisson models via 
quasi-likelihood inference [36, 4]. We used the software DMU version 6.0, 
release 5.1 [23, 24], which implements multivariate generalized linear mixed 
models efficiently and facilities to specify the required covariance structure 
of the random components. 

3.2 A Poisson approximation for discrete- time models 

The inference for the discrete time models presented here is equivalent to 
making inference for a Bernoulli model where the response variable is a vari- 
able indicating whether a culling event occurred. Here, we present a Poisson 
approximation that will be useful for performing likelihood-based inference 
for the discrete time models, specially in the case of very large populations 
with relatively low occurrence of culling. More precisely, for each time t 
(t = 0,1. ... ,t, where r is the maximum observed time) and for each in- 
dividual i (i = 0, 1, . . . ,n t where n t is the number of individuals at risk in 
the time t), we constructed a pseudo observation of a response variable Y it 
taking the value 1 if the individual % died at time t and otherwise. Ex- 
ploring a coincidence of the likelihood function, we treated the observations 

Y it , for t — 0, 1 , r and i — 0, 1, . . . , n t , as independent Bernoulli random 

variables with P(Y it = 1) = p it , where p it is the hazard probability for the 
individual i at the time t. 

The following condition will ensure that the likelihood function of a suit- 
ably defined Poisson model will approximate the likelihood function of the 
Bernoulli models considered here. Suppose that 

max p t , nt -> , (3) 

1<I<T 

where p t>nt = maxi<;< nt p it , and 

nt 

Pu = 7* (fixed) for all t as min n t — > oo . (4) 
The general result on convergence for arrays of probabilities presented 
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in [7] implies that, under the conditions (3) and (4), for each time t and 
r = Q > 1 >---> 

lim P(W nt = r) = ^f- , (5) 

n t —>-0 r\ 

where W nt = X^=i ^*t- Note that the right side of the equation (5) is the 
density probability function of a Poisson random variable with expected value 
equal to 7$. In other words, if the number of individuals at risk in each time t 
(i.e. n t ) is larger enough and the probability of death in each time t is small 
enough, then the likelihood function of the Bernoulli model is approximated 
by the likelihood function of a Poisson model. 

Note that the conditions given in (3) and (4) are reasonable in the ex- 
ample presented here, since the number of individuals is of the order of 100 
thousands at the first observed period and is of the order of 16 thousands 
at the last observed period. Moreover, for a given individual i, the pseudo 
random variables Y it (t — 1, . . . , tj) is a sequence of zeros until time — 1 
and takes the valuel at £j only if this individual was culled, therefore p it will 
be small in general. 

3.3 Continuous models with stratification 

The examples of survival models with continuous time presented here re- 
quired the use of stratification of the baseline function, which is a standard 
technique of survival analysis. For example, we used the number of days 
from the first parity to the culling day as one of the characterizations of the 
longevity of animals (sows or dairy cows), in the two examples presented 
here. Exploratory analyses of the data used indicated that there were differ- 
ences in the mortality rates between the different parities and that the form 
of the baseline functions were not the same for the different parities, i. e. the 
baseline functions for the different parities were not proportional. A way to 
circumvent this issue is to use the stratification technique defined here. This 
type of model is referred in the literature of dairy cattle longevity as the 
"lactation basis model" (see [29]). 

Define t Pi as the time of the occurrence of the p th parity of the individual 
i (i = 1, . . . , n and p — 1, . . . , Pi where Pi is the maximum parity observed 
for the individual i) and t* as the observed time (death or censure) for the i th 
individual. The hazard function for the i th individual at the p th parity, for 
% = 1,2, ... , n and for p = 1, . . . , Pj conditional on the random components 
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U = u and V = v is then 

\ p {t\U, V) = Y i!P (t)\ p (t) exp (X itP (3 + Z iU + W,v) , 

for t G [t Pi , tr Pi+ i)) when p < Pi and for t G [tp v t*) when p = Pi. Here A p (-) 
is the baseline hazard function for the p th parity. The variable Yi :P (t) is equal 
to 1 if the individual was at risk at t G [t Pi ,t( Pi +i)) and otherwise. In this 
study, the time dependent explanatory variables will be a function of p and 
then X ijP (t) = X i:P for all t > 0. 

3.4 Comparison of the statistical information between 
discrete and continuous time models 

When applying the longevity models described above it is often possible to 
choose between models defined with continuous time and models defined 
with discrete time. For instance, in the examples studied in section 5 the 
longevity can be characterized by the time between the first parity and the 
culling (continuous time) using a stratification by parity or, alternatively, 
by the number of survived parities (discrete time). The complexity and the 
statistical stability of these two types of models differ, which is an important 
aspect to take into account when deciding which type of model to use. We will 
argue next that the likelihood function of a continuous time model is typically 
much flatter than the likelihood function of a corresponding discrete time 
model. This explains certain anomalies in the inference under continuous 
time models. 

In order to simplify the exposition, we assume a scenario where two mod- 
els are considered: a discrete time model where the number of survived par- 
ities is modeled and a continuous time model where the number of survived 
days is modeled but there is a stratification by parity (as described in sec- 
tion 3.3). Moreover, we restrict the discussion to the case where there is 
only one cause of culling and there is only one random component, say U, 
in the model. Here, U will be a multivariate normal random variable with 
E(XJ) = 0, Var(XJ) = Acr 2 , where A is a known matrix. We denote the 
density function of the distribution of U by $ (•, a 2 ) and represent the mul- 
tiple integral for integration with respect to this density function by a single 
integration sign. 

The likelihood function of the piece- wise constant hazard model stratified 
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by parity based on the number of survival days is 



J Yl exp {-Aipfc exp {riipk)} {A ipk exp (rj ipk )} yipk $ (u, a 2 ) du , (6) 

i,p,k 

where the indices i,p and k indexes the individual, the parity and the time 
intervals in each periods, respectively. Here, yi P k is a indicator variable taking 
the value 1 if the i th individual died at the p th parity in the k th time interval 
and otherwise; A ipk is the length of the k th time interval at the p th parity; 
and rji P k = \og(\ p k) + X ipfc /3 + ZjU is the linear predictor. 

On the other hand, by the Poisson approximation discussed in section 
3.3, the likelihood function of the discrete relative risk model based on the 
number of survived parities is approximated by 

/ Y[ exp {- exp (m P )} {exp {r] vp )} y ^ $ (u, a 2 ) du , (7) 

i,p 

where % and p indexes the individual and the parity, respectively; yi P is the 
indicator variable that the i th individual died at the p th parity and r\ ip = 
log(Ap) + X ip /3 + Z^u is the linear predictor. 

Note that, when there are many censured observations and when there 
are many survived parities or many observed time intervals in each parity, the 
variables y ip k and y ip will take the value zero in most of the cases. Therefore, 
under these conditions, the integrand in (6) will be dominated by 



JJ exp {-A ipk exp {i] ip k)} = exp <^ - ^ A ipk exp (i] ipk ) > 

i,p,k I i,p,k ) 



(8) 



and the integrand in (7) will be dominated by 

n exp {- exp (r} ip )} = exp <^ - exp (r} ip ) \ . (9) 



%,p 



Note that the sum in the right side of (8) has typically a larger number 
of parcels as compared to the sum in the right side of (9). Moreover, the 
parcels in the sum in the right side of (8) are multiplied by the factors 
A ip k which are typically larger than 1. Note, moreover, that if one reduces 



14 



the size of the intervals used in the piecewise constant hazard model, the 
factors Ajpfc decrease, but the number of parcels in the sum in the right 
side of (8) increases; on the other hand, if one decreases the number of 
intervals, the number of parcels decreases, but the factors A ipk might increase 
substantially. This makes the likelihood function of the piecewise constant 
hazard models much flatter, since the likelihood function is positive. More 
precisely, the curvature of the support curve (i.e. the graph of the logarithm 
of the likelihood function) near the maximum likelihood estimate will be 
much smaller for the piecewise constant hazard models. Therefore, the Fisher 
information will be smaller, and the maximization of the likelihood function 
will be a much harder problem, for the piecewise constant hazard models as 
compared to the corresponding discrete relative risk model. 



4 Decomposition of the phenotypic variance 
and heritability 

The notion of heritability is crucial in quantitative genetic applications; it 
measures the magnitude of the detected additive genetic signal relative to 
the total variance of a trait of interest and is typically used to quantify the 
potential response to selection. Heritability (in the narrow sense) is defined 
operationally by the ratio 



where a 1 is the variance of an unobserved random component representing 
all the additive genetic variation (termed the genetic variance) and a 2 is the 
variance of the trait of interest (called the phenotypic variance). For details 
see [22, 18]. In the classic scenario, where gaussian linear mixed models are 
used, the calculation of the heritability is straightforward because the pheno- 
typic variance can be expressed as the sum of the genetic variance a 2 g and the 
variances of the other random components present in the model. However, 
the calculation of the heritability is much more complicated in the context of 
survival analysis considered here and are, therefore, developed in full details 
below. It will be necessary to introduce some details of the counting process 
behind the models used and to define the trait of interest carefully. Here, 
we present the main results informally and expose the precise formal defini- 
tions and the technical details using the mathematical machinery of counting 
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processes in appendix A. For simplicity of exposition, we discuss only the 
one-dimensional case were only one cause of culling is present (referred as 
death) and assume the same model structure presented in section 3 but with 
dimension one. 

Let T be a non-negative random variable representing the observed sur- 
vival time and D be an indicator variable taking the value 1 if a death 
is observed and otherwise (censure). The random variable T will take 
values in R + if the time is continuous or in Z + {i.e. in {0,1,2,...}) if 
the time is discrete. We assume that for each time t a vector X(t) of 
explanatory variables (possibly changing with time) is known. Denote by 
X(£) = {X(s) : < s < t} the trajectory (or path) function until the time 
just before t (denoted by t—). Additionally, we consider two independent 
gaussian random components, U and V, representing the additive genetic 
effects and some environmental effects, respectively. Here, we follow the 
same convention as in section 2, so Var(U) = Acr| and Var(V) = Iof, 
where A is the relationship matrix and I is an identity matrix with suitable 
dimensions. The hazard function for the i th individual takes then the form 

A,(t|U,V) = A (t)exp(x;(t)/3 + Z;u + W;v) . (11) 

In order to describe the models presented in section 3 properly, we need 
to introduce two stochastic processes: the death (or event) counting process 
and the at risk process. To define the death process we consider, for each 
t, the random variable N(t), which indicates whether a death occurred until 
the time t, i.e. N(t) = 1(T < t, D = 1). The stochastic process N = 
{N(t) : t G R+ or Z + } is the death counting process. It is convenient to 
introduce the notation dN(t) to denote the increments of the counting process 
N at each time t, i.e. dN(t) = N(t) - N(t-), where N(t-) = lim A +o N(t - 
A). Furthermore, the at risk process given by Y = {Y(t) : t G K+ or Z + }, 
where, for each time t, the random variable Y(t) takes the value 1 when the 
individual is at risk at time t or otherwise. 

We will calculate the phenotypic variance and the heritability for two 
characteristics of interest: the hazard function evaluated at a given time and 
the accumulated hazard function up to a given time. These two charac- 
teristics are associated to the intensity and the cumulative intensity of the 
counting process of death N. Due to the structure of the models used here 
(see equation (1) and (11) for instance), the entire hazard function of an 
individual is multiplied by the factor exp(u), where u is the realization of 
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the genetic random component U for that individual; so, according to this 
model, when selecting individuals by the predicted values of U (as it is the 
usual practice in animal breeding), the entire hazard curve is affected and 
therefore, it is indifferent if one measures the potential response to selection 
(i.e. the heritability) of the hazard evaluated at a given time-point, say t, or 
even of the accumulated hazard up to t. 

In order to study the heritability related to the hazard function evaluated 
at a given time t (in M + or in Z+), we define the conditional random variable 

y(t) = dN(t)\T>t,X(t). (12) 

Note that E (y(t) | U = u, V = v) = A (t | U = u, V = v), therefore y(t) 
will be treated as the trait of interest, for which we calculate the phenotypic 
variance and the heritability. A more precise definition of y(t) using counting 
processes theory is given in appendix A. There we show, taking advantage of 
the basic theory of martingales and using a suitable Taylor expansion, that 
in the continuous time case the variance of y(t) is approximated by 

Var [y(t)\ « Y(t) {[X*(t)f (a 2 g + * 2 e ) + 0A*(i)} , (13) 

where X*(t) is the hazard function evaluated at t with vanishing random 
components, i.e. 

A*(t) = A(t|U = 0,V = 0) = A (t)exp{x;(t)/3} . (14) 

Note that according to (13) the variance of y(t) is additively decomposed 
in three components: one depending on a g , one depending on a\ and one 
depending on 0. This is a situation analogous to the classic gaussian linear 
mixed models. Here, the component of the variance of y(t) associated to a g , 
namely [Y(t)X*(t)} 2 a g , plays the role of the genetic variance and therefore 
the heritability for the hazard evaluated at the time t is given by 

%m*nt) - "} , ■ (is) 

9 e A*(t) 

Analogous calculations (see appendix A) yield for the discrete time case 
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that the variance of y(t) is approximated by 



Var [y(t)\ « Y(t) {[A*(i)] 2 {a 2 g + a 2 e ) + <f>\*(t)[l - X*(t)]} . (16) 

Here, as in the case of the continuous time, the variance of y(t) is additively 
decomposed in three separate a components depending on a 2 , on a 2 and on 
0, respectively. Therefore, the heritability is approximated by 

2 

hit) ~ Y(t) ^ . (17) 

°S t0 e + ^1-A*(t) 



Moreover, when a Poisson approximation is used, the approximated heri- 
tability takes the form given in (15). 

We turn now to the calculation of the heritability of the accumulated 
hazard function evaluated at a time point t. In the continuous time case we 
define the accumulated hazard function, for each t G R+, by 

A(t) = [ Y(s)X(s)ds, (18) 
Jo 

and in the discrete time case the accumulated hazard function is given, for 
each t £ Z + , by 

A(t) = £V(«)A(*) , (19) 

s<t 

where A(-) is the hazard function in play. Furthermore, the conditional ran- 
dom variable 

Z(t) = N(t)\T >t,X(t) (20) 

is such that E (Z(t) \ U = u, V = v) = A (t \ U = u, V = v). Therefore, 
Z(t) will be treated as the trait of interest, for which we calculate the phe- 
notypic variance and the heritability. A more precise definition of Z(t) using 
counting processes theory is given in appendix A. It can be shown (see 
appendix A) that in the case of the continuous time the variance of Z(-) 
evaluated at a given t G K+ is approximated by 

Var [Z(t)\ » [A*(t)f {a 2 g + a 2 e ) + 0A*(i) (21) 
where A*(t) = f*Y(s)\*(s)ds. Therefore, the heritability is approximated 
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by 



/ '' ,/, ^ + ^ + 0[A*(*)]" 1 • (22) 



In the case of the discrete time the variance of Z(-) evaluated at a given 
t G Z + is approximated by 

Var [Z(t)\ « [A*(t)] 2 (a 2 + a e 2 ) + 7 (t), (23) 

where 7(4) = ^* =0 <pY {t)X* {t)[l - \*{t)\ and the cumulative hazard is A*(t) = 
^* =0 F(s)A*(s) .Then, the heritability is approximated by 

^<*> * ff 2 + ^ + ^y(*)[A*(*)]- 2 ' (24) 

When the Poisson approximation is used the heritability is approximated by 
the right side of (22). 

The calculations above are performed for one individual. We propose 
that given a population of n individuals, the phenotypic variances and the 
heritabilities for the hazard function or the cumulative hazard (evaluated at 
a time t) should be calculated using the estimated variance components and 
an estimated hazard or an estimated cumulative hazard as described below. 
In the continuous time case, where the PCHM is used, the observed time is 
split in intervals I k = [t k , ifc+i), for k — 0, . . . , K where — t < t± < . . . < 
tx < tx+i = t, where r is the largest observed time. Given that U = u 
and V = v, where u and v are the BLUP of the random components, the 
conditional hazard is defined, for each t G Ik, by 

x*(t) = \ k exp [x; >fc /3 + z; u + w; v ] = exp [ Vhk (t)] , 

where r] ijk (t) = log( X k ) + X ' ik (3 + Z -u + W -v and X ijfe = Xj(i) for t G I k . We 

obtain f]^ k (t) = log(Afc) + X, k f5 + Z^u + W^v from the adjusted model and 
then calculate, for each time t G I k , 

\*(t) = X* k = exp [fj k ] , 
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where r] k = n 1 Yli=i Vi,k- The cumulative hazard is estimated, for t £ by 



k-i 



3=0 

where Aj = tj + i — tj is the length of the interval Ij. 

In the discrete time case the hazard, conditionally that U = u and V = v, 
is given, for each t £ Z + , by 

X*(t) = A t exp [x;(0/3 + Z; U + Wjvl = exp [ Vi (t)] , 



where ^(t) = log(A t ) + X^(t)y9 + Z^u + W-v. We obtain ^(t) = log(A t ) + 
X-(t)/3 + Z-u + W-v from the adjusted model and then calculate for each 
timet (f = 0,1,2,...,) 

~X(t) = exp [fj(t)} , 

where fj(t) = n^ 1 YH=i Vi{t)- The cumulative hazard is estimated by A*(t) = 
Y^j=o ^*0')- The heritabilities will be estimated at X*(t m ) and A*(t m ) where 



t m represents the median survival time, i.e. P{T > t r 
tive examples. 



0.5 in the illustra- 



5 Two illustrative examples 

The methods presented above will be illustrated by two examples: one in- 
volving the longevity of sows (section 5.1) and another characterizing the 
longevity of dairy cattle (section 5.2). In both examples, we will characterize 
the longevity of female animals in two different ways: first, by the length of 
the productive life of the animals (sows or cows) measured by the number of 
days from the first parity to the culling day (ND); second, by the number of 
survived parities (NP). Models based on continuous time and models based 
on discrete time were applied to study ND and NP, respectively. 

Incomplete observations are present in both examples because some ani- 
mals were still alive at the end of the study, some animals were moved, sold 
or exported during the study and due to intentional right truncation. In the 
study of sow longevity the proportion of censure was relatively low (9.6%), 
while in the study of longevity of cows a much higher proportion of censure 
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was observed (23.5%). Furthermore, in the study of longevity of cows it 
was possible to distinguish two culling causes: slaughter (67.7% of the ob- 
servations) and death (only 8.8% of the observations), which characterizes a 
typical scenario of competing risks. Therefore, it is natural to use bivariate 
models in this example to model these two culling causes simultaneously. 
Bivariate models will be used in the example of sows longevity to model ND 
and NP simultaneously and use that to discuss how much the information of 
these two characterizations of the longevity overlap. 

5.1 Longevity of sows 

The data in this example was retrieved from the registers of the Danish 
Pig Center of 200,084 pure Landrace sows that were giving birth in the 
years between 1999 and 2010. Only the sows that had the first parity in 
the comprised period were included in the study. All the models proposed 
for ND and NP included the following explanatory variables: age at the 
first parity, year and season at the first parity, litter size (total number of 
piglets born per parity) and herd year size (total number of sows farrowing 
per herd per year). The litter size and the herd year size were both time 
dependent variables. In addition two random components were included: a 
sire component with pedigree, representing the sire additive genetic effects 
and a herd-year component representing the environment effects. Since a sire 
model was used, the additive genetic variance was estimated by multiplying 
the additive sire variance by four. 

The Kaplan-Meyer estimate of the median survived NP was 2 parities 
(CT95% 2;3) and of the median survived ND was 320 days (95% CI, 319 
- 322). The median number of days between two consecutive parities for 
the sows that did not die in the respective period varied from 152 to 155 
days. Thus, the median survived time of 320 days from the first parity is 
approximately the time necessary for the sow, give it had one parity, have 
the second and survive until a third parity. 

Figure 1 displays the logarithm of the cumulative hazard curves for ND 
stratified by parity. The curves per parity presented similar behavior, how- 
ever, they are not parallels indicating a non proportional effect of parity. 
The dashed lines represents the cut-points of the intervals where the baseline 
hazard function was assumed to be constant. 

Table 1 presents the estimates of the variance components and heritabil- 
ities from the piece-wise constant hazard model stratified by parity based 
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on the ND and for two discrete relative risk models based on the NP (an 
exact model and an approximated model via Poisson approximation). The 
estimates of variance for the sire and the herd-year components were very 
similar among the three adjusted models (~ 0.065 with SE ~ 0.003 for the 
sire component and ~ 0.220 with SE ~ 0.009 for the herd- year component). 
Comparing the execution time between the two adjusted discrete models the 
approximated models was around 72 times faster compared to the time spent 
to execute the exact model. 

A bivariate model describing the ND and the NP traits (using a Poisson 
approximation for the discrete marginal models) was fitted (see Table 2). 

Table 1: Estimated variance components (with asymptotic standard error 
in parenthesis) and heritabilities for the number of days (ND) and for the 
number of parities (NP) at the sows longevity study. 



Source 


ND-PCHM a 


NP-DRRM b 


NP-DRRMp 


Sire 


0.071 (0.005) 


0.066 (0.002) 


0.063 (0.003) 


Herd- Year 


0.217 (0.010) 


0.220 (0.005) 


0.219 (0.009) 


Dispersion 


3.089 (0.002) 


0.978 (0.002) 


0.653 (0.001) 


hi* 




0.098 


0.096 


hi* 


0.180 


0.161 


0.165 


Execution time c 


48.5 min 


8h 28 min 


7 min 



a PCHM: Piece-wise constant hazard model stratified by parity. 
b DRRM: Discrete relative risk model. 

c DRRMp: Discrete relative risk model via Poisson approximation. 
d Estimated heritability for the hazard (A) at the median survival time. 
c Estimated heritability for the cumulative hazard (A) at the median 
survival time. 

f Execution time (Intel i5 processor). 
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Table 2: Estimated variance components and correlations from the bivariate 
model based on the number of days (ND) and number of parity (NP) at the 



sows longevity 


study. 






Source 


Trait 


a 2 (SE) 


Cor (SE) 


Sire 


ND 
NP 


0.022 (0.001) 
0.081 (0.004) 


1.000 (0.004) 


Herd- Year 


ND 
NP 


0.332 (0.015) 
0.231 (0.010) 


1.000 (0.002) 


Dispersion 


ND 
NP 


3.130 (0.002) 
0.657 (0.001) 








1 - - 2 ••• 


Parity 

3 - 4 5 •-• 



1 -I 




1 2 5 10 20 50 100 



Figure 1: Log of cumulative hazard curves for the number of days trait (ND) 
stratified by parity at the sow's longevity study. Axis x is presented on the 
logarithm of days scale. The vertical dashed lines are the cut points for the 
intervals where the baseline function was assumed to be constant. 
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5.2 Longevity of dairy cattle 

The data analyzed in this example was retrieved from the register of Knowl- 
edge Centre for Agriculture of 142,133 Jersey Danish dairy cows that were 
calving in the period between 1990 to 2006. The explanatory variables (fixed 
effects) included in the models were: herd year size (number of calving per 
herd per year), year and season of each parity, as time dependent variables 
and age at the first parity and a coefficient of general heterosis (defined as 
a continuous explanatory variables), as time independent variables. Two 
random components were included in the model, a combination of herd- year- 
season component representing the environment effects and a sire component 
with the pedigree representing the genetic effect. In addition, the cows could 
be culled by one of two possible reasons: death and slaughter. 

The Kaplan-Meyer estimate of the overall median NP was 3 parities (95% 
CI, 3;3) and overall median ND was 826 days (95% CI, 820 - 831). The 
median number of days between 2 consecutive parities for the cows that did 
not die in the respective period varied from 361 days to 367 days. Figure 
2 displays the stratified cumulative incidence probability curves for ND for 
the slaughter and the death, respectively. The dashed lines represent the 
cut points for the intervals where the baseline hazard were assumed to be 
constant for both specific reasons. 

Tables 3 displays the estimates of the variance and covariance components 
from the competing risk models describing the ND and NP traits. The sire 
end herd-year-season variances for death rate based on the ND trait were very 
small, both for the death and the slaughter rate. In contrast, the discrete time 
model detected significantly larger variances of the random components. The 
failure of the continuous time models can be explained by the large amount 
of censure in the data. 
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Table 3: Estimated variances and covariances components for the number 
of parities (NP) and the number of days (ND) from a discrete and continuous 
multivariate competing risk model, respectively, at the dairy cattle longevity 
study. 



Source 


Cause 


NP 

o 2 (SE) 


Cov (SE) 


ND 

a 2 (SE) 


Cov (SE) 


Sire 


Death 
Slaughter 


0.102 (0.041) 
0.029 (0.003) 


0.002 
(0.004) 


0.7xl0" 3 (0.006) 
0.008 (0.002) 


0.002 
(0.003) 


HYS a 


Death 
Slaughter 


0.458 (0.015) 
0.061 (0.002) 


-0.055 
(0.004) 


0.1 xl0~ 5 (0.079) 
0.3xl0- 6 (0.010) 


0.7xl0~ 7 
(0.021) 




Death 
Slaughter 


0.691 (0.002) 
0.710 (0.002) 




5.165 (0.005) 
4.148 (0.004) 




hi- 


Death 
Slaughter 


0.013 
0.040 










Death 
Slaughter 


0.044 
0.108 




0.2x10" 
0.012 


3 



a Herd-year-season random component. 
b Marginal dispersion parameters. 

c Marginal heritabilities for the hazard (A) at the median overall survival 
time. 

d Marginal heritabilities for the cumulative hazard (A) at the median 
overall survival time. 
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Slaughter 



Death 




Figure 2: Cumulative incidence probability curves for the ND stratified by 
parity for the cow's longevity study. The dashed vertical lines are cut points 
for the intervals were the hazard was assumed constant. Pj(t) = P(T < 
t, J = j) for J = (Death, Slaughter). 
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6 Discussion 



The several variants of models for studying longevity in quantitative genetics 
can be described using the following general approach. The models are one- 
dimensional and describe the conditional hazard function, given the random 
components, as taking the form 

\(t\V =u,V = v) = A (t) exp (x;(0/3 + Z[u + Wjv) , 

for t in R + or in Z + . The models based on continuous time used in the 
literature differ essentially in the form of the baseline hazard function Ao(-) : 
the classical Cox model (with frailties) assumes the hazard function to be 
a smooth function (almost everywhere) imposing in this way almost no re- 
striction on A (-) (see [19, 26, 28]); the piece-wise Weibull models (see [9]) 
assume A (-) to be a saw function (i.e. a piecewise linear function); while 
the piece-wise constant hazard models (as described here) assume Ao(-) to 
be a step function (i.e. a piece-wise constant function). These models are 
nested, in the sense that saw functions are smooth almost everywhere, and 
step functions are particular cases of saw functions. Moreover, smooth func- 
tions can be arbitrarily approximated (point-wisely) by saw functions or by 
step functions. A fourth category of models, which was introduced here, are 
the piece-wise constant hazard models with free dispersion parameter. These 
models are an instance of quasi-likelihood based models (see [36, 4]). They 
contain the traditional piece- wise constant hazard models as a particular case 
(by setting the dispersion parameter to be constant equal to 1) and might 
arbitrarily approximate the piece-wise Weibull models and the Cox frailty 
models. Although in many practical situations those models yield equivalent 
representations of the data and produce similar results, this is not always the 
case for the quantitative genetics applications discussed here. 

One of the interesting features of the piece-wise constant hazard model 
with (free) dispersion parameter presented here is that the phenotypic vari- 
ance can be approximately decomposed as the sum of parcels representing 
the additive genetic effects, environmental effects and unspecified sources of 
variability. This is analogous to the classic decomposition of the phenotypic 
variance obtained in gaussian linear mixed genetic models. Here, the compo- 
nent of the phenotypic variance involving the dispersion parameter plays 
the rule of the residual variance in the classic models. Note that, if we fix 
the dispersion parameter 0, say by setting it equal to 1 as in the piece-wise 
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constant hazard model, it might be that we obtain models that describe the 
data well and the referred decomposition phenotypic variance would still ex- 
ists, but this would be equivalent to set the residual variance to be constant, 
which clearly causes problems in the representation of quantitative genetic 
phenomena. For instance, in a sire model the genetic variance transmitted 
by the dam is not captured by the random component representing the sire, 
therefore, this part of the variability should be represented in the residual 
variance (i.e. should compose part of the residual variance); but this is not 
possible to occur if the residual variance is set to be constant. This issue oc- 
curs also in the piece-wise Weibull models (see [9]). Furthermore, in the most 
general case where A is only assumed to be smooth (the Cox frailty model) 
it might be argued that Ao could be point wisely approximated (almost every 
where) by a sequence of step functions obtained even from distributions in a 
model with dispersion parameter, which would have a proper decomposition 
of the phenotypic variance. However, in this case it would not be possible to 
identify the dispersion parameter (as essentially argued in [13]); therefore, 
one would still have the representation issues referred above. 

The piece-wise constant hazard model and the piece-wise Weibull model 
are parametric models, making it possible to implement relatively efficient 
inference techniques. This explains why they are implementable for complex 
and large models as in quantitative genetics typical applications. Moreover, 
when the baseline hazard function A (-) is assumed to be piece-wise linear (a 
saw function), as in the piece- wise Weilbull models, the rate of approximation 
of Ao(-) to a smooth function (as in the frailty model) is improved (specially 
if there are regions where the smooth baseline is steep) as well argued in 
[9]; however, in this case the inference via generalized linear mixed models 
would not be possible without using profile likelihood techniques (to estimate 
the shape parameter), which rules out the direct use of extensions contain- 
ing a dispersion parameter or multivariate extensions already implemented. 
Finally, a disadvantage of the piece-wise constant hazard models and the 
piece-wise Weibull models is that both depend on an arbitrary choice of time 
cutting points. 

Models based on discrete time (i.e. variants of the discrete relative risk 
models) has been occasionally used in quantitative genetic previously; how- 
ever, without incorporating a dispersion parameter. These models do not 
assume any pre-specified form for the baseline hazard function and also not 
require any arbitrary definition of time cutting points; therefore, they are 
similar to the Cox frailty models. Due to the coincidence with the Bernoulli 
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models with random components (with or without dispersion parameter), 
techniques for parameter inference are at hand, allowing for efficient imple- 
mentations capable to handle complex large problems as illustrated here. 
Moreover, these models allow to incorporate time- dependent variables (by 
splitting data with record for each survived time for each individual) both 
as fixed effects and as random components. The Poisson approximation pre- 
sented here might represent a substantial save in computational resources 
and yields essentially the same results as the exact inference, as illustrated 
in the example of sows longevity presented in section 5.1. 

The continuous time models and the discrete time models yielded simi- 
lar results in the example of sows longevity. However, the continuous time 
models failed in the example of dairy cattle longevity (section 5.2), which 
can be explained by the fact that, in this case, both rare events and high 
frequencies of censored observations are observed. This is in line with the 
limitations of the continuous time models discussed in section 3.4 and shows 
that those limitations do occur in practical situations. 

The lack of multivariate techniques for analyzing several traits simulta- 
neously has been a serious limitation in the application of survival models in 
quantitative genetics [30, 11]. This lack not only makes the study of the rela- 
tion between longevity and other traits impossible, but also limits very much 
the possibility of performing time-to-event-analysis in the presence of com- 
peting risks. Indeed, if there are two competing culling reasons, one might 
use standard survival models for modeling the rate of occurrence of one of the 
culling reasons by considering the observations where the other culling reason 
occurred as censored. However, this naive approach implies in the implicit 
use of the assumption that the competing culling reasons are independent. 
The example of dairy cattle longevity presented here illustrates that this as- 
sumption might be not reasonable. Indeed, in this example we detected a 
significant negative correlation between the herd-year-season random com- 
ponents; moreover, a range of fixed effects were found to be significant for 
both culling reasons, implying that a naive use of marginal models would 
violate the assumption of absence of informative censoring. 
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A Technical details on the decomposition of 
the phenotypic variance 

Let T be a non-negative random variable representing the observed survival 
time and D be an indicator variable taking the value 1 if a death is observed 
and otherwise (censure). Denote by X(i) = {X(s) : < s < t} the tra- 
jectory function of the known vector of explanatory variables until the time 
just before t. Additionally, we consider two independent gaussian random 
components, U and V, representing the additive genetic effects and some 
environmental effects, respectively. 

We define the death counting process N = {N(t) : t G R+ or Z + }, where, 
for each time t, N(t) = 1(T < t,D — 1). The increment of the death 
process at the time t is given by dN(t) = N(t) — N(t—), where N(t—) = 
lim^o N(t — A). The at risk process is given by Y = {Y(t) : t G M+ or Z + }, 
where, for each time t, the random variable Y(t) takes the value 1 when the 
individual is at risk at time t or otherwise. 

The death counting process and the at risk process are both adapted 
to the filtration F = {J~t- '■ t G R+ or Z + }, where, for each time t, Tt- = 
a {N(s), Y(s),X(s), < s < t} is the a-algebra generated by N(-), Y(-), X(-) 
up to t— (i.e. up to the time just before t) . Here, T t - represents what is 
known up to (but not including) time t. 

Conditional on the random components U and V, under the standard 
independent censoring assumption (see [3, 19]), we have that 

P[dN(t) = 1|U = u,V = v,7i-] = F(t)A(t|u,v), (25) 

where A(f|u,v) = X (t) exp [X(t)'/3 + Zu + Ww] and Y(t)\(t\u,v) is the 
intensity associated to the counting process N at time t. Moreover, the 
cumulative intensity of N up to time t is A(t|u, v) = J ' F(s)A(s|u, v)ds. 
Note that A = {A(t|u,v) : t G R + or Z + } is predictable with respect to the 
filtration F . Define, for each time t, 

M(t) = N(t) - A(t) . (26) 

The process M = {M(t) : t G R+ or Z + } is a martingale. Here, (26) is 
the Doob-Meyer decomposition of the death counting process where A is the 
compensator. Define, for each time t, the increment of the martingale M by 
dM(t) = M(t) - M(t-), where M(f-) = lim A40 M(t - A). Since M is a 
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martingale, E[dM (t)\J-t-\ = for each time t, which shows that M plays the 
role of the residuals of the model. The predictable variation process of the 
martingale M at the time t is 

(M) (t) = [ Var[dM(s)\F s -\ 
Jo 

and 

d(M) (t) = Var[dM(t)\F t -\ . 

Now, we study the decomposition of the phenotypic variance related to 
the hazard function evaluated at a given time t. Define the conditional 
random variable 

y(t) = dN(t)\T t - , (27) 

which is equivalent to y(t) informally defined by (12) in section 4. 
Given the random components U and V, 

E[y(t)\V,V\=Y(t)\(t\u,v) (28) 

and 

Var [y{t) |U, V] = (f)Var [dM(t)\XJ, V, T t -\ , (29) 

where <fi is the dispersion parameter. Using (28) and (29), the variance of 
y(t) is given by 

Var[y(t)\ = Var {E [y(t) |U, V] } + E {Var [y(t)\U,V]} (30) 
= Var{Y(t)X(t\V,V)} + E{(j)Var[dM(t)\V,V,T t -]} . 

A Taylor expansion argument (see [16], page 118) yields 

E{Var [y(t)\U,V]} w Var [Y(t) \ U = 0, V = 0] 

= 4>Var[dM(t)\XJ = 0,V = 0,^-] . (31) 

Moreover, a second order Taylor expansion (see [16], page 118) yields 

Var{E[y(t)\XJ,V}} « o\ \^-E [y(t)\V, V\Y 

I du J (U=o,V=o) 
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+ CT e{|-^W)|U,V]j 2 

I dv J (U=o,V=o) 

= [Y(t)\*(t)] 2 (al + a 2 e ) , (32) 

where u and v are realisations of the random components. Moreover, A* (t) = 
X(t\V = 0,V = 0). 

In the continuous time case, for each t G R+, 

Var[dM(t)\U = 0,V = 0, F t -] =Y(t)\*(t) , 

(see [19, 3]) and therefore, using (31) and (32) yields 

Var [y(t)\ « Y(t) {[A*(t)] 2 (a 2 g + a 2 ) + 4>\*{t)} , 

which proves the decomposition of the phenotypic variance given in (13). 
In the discrete time case, for each t G Z + , 

^ar[dM(£)|U = 0,V = 0, T t -\ = Y(t)X*(t) [1 - X*(t)\ , 

(see [19, 3]) and therefore, using (31) and (32) yields 

Var [y(t)} « F(0 {[A*(t)] 2 + a e 2 ) + 0A*(t) [1 - A*(0]} , (33) 

which proves the decomposition of the phenotypic variance given in (16). 

For the decomposition of the phenotypic variance related to the cumu- 
lative hazard function evaluated at a given time t we define the conditional 
random variable 

Z(t) = N(t)\T t - . 

This is equivalent to conditional random variable Z(t) informally defined in 
(20). Note that, 

E [Z{t)\\J, V] = A(t|U, V) and Var [Z{t)\\J, V] = (M) (t) |U, V . 

The decomposition of the variance of Z(t) is obtained by applying an analo- 
gous calculation as in the decomposition of the variance of y{t). This yields 
the equations (21) and (23) for the decomposition of the phenotypic variance 
for the continuous time and the discrete time cases. 
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